Local Elastic Constants in Thin Films of an FCC Crystal 
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In this work we present a formalism for the calculation of the local elastic constants in inhomoge- 
neous systems based on a method of planes. Unlike previous work, this formalism does not require 
the partitioning of the system into a set of finite volumes over which average elastic constants are 
calculated. Results for the calculation of the local elastic constants of a nearest neighbor Lennard- 
Jones fee crystal in the bulk and in a thin film are presented. The local constants are calculated at 
exact planes of the (001) face of the crystal. The average elastic constants of the bulk system are 
also computed and are consistent with the local constants. Additionally we present the local stress 
profiles in the thin film when a small uniaxial strain is applied. The resulting stress profile compares 
favorably with the stress profile predicted via the local elastic constants. The surface melting of a 
model for argon for which experimental and simulation data are available is also studied within the 
framework of this formalism. 

PACS numbers: 68.08.-p,68.35.Gy,68.60.-p 
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I. INTRODUCTION 

Knowledge of the local elastic constants in inhomo- 
geneous systems is of significant theoretical, experimen- 
tal, and industrial interest. As nanofabrication technolo- 
gies improve and allow for the design and construction 
of nanoscopic devices, understanding the mechanical re- 
sponse of materials at nanometer length scales will be- 
come increasingly important. In particular, deviations 
from bulk, continuum behavior may lead to complica- 
tions in the manufacturing of such devices. For example, 
in the microelectronics industry, the mechanical collapse 
of photoresist structures below lOOnm may limit the ul- 
timate density of memory storage devices or the perfor- 
mance of microprocessors H, |^ . 

In nanoscopic structures, interfaces are likely to play a 
major role in apparent deviations from bulk, continuum 
behavior Q. The interface could either weaken or re- 
inforce the overall mechanical behavior of the structure, 
depending on the nature of interactions between adjacent 
domains and the size of the structure. Understanding 
how mechanical properties vary near interfaces or free 
surfaces would provide insights into such phenomena. 

Knowledge of interfacial behavior is crucial for under- 
standing the adhesion of thin polymer films, where the 
inter-diffusion of the polymers and the molecular mobil- 
ity near the film boundaries play a significant role 
Properties such as adhesion, dewetting, and surface melt- 
ing in thin films are likely to be controlled by processes 
that occur within the first few nanometers of the inter- 
face. It would therefore be beneficial to have the ability 
to measure (computationally or experimentally) physical 
properties with molecular spatial resolution. 

A microscopic definition for local elastic constants has 
been proposed in the literature Q . Implementation of 
that formalism requires that layer- averaged local elastic 
constants be determined. For inhomogeneous systems, 
the results from averaging over a particular layer depend 
strongly on the size and position of that layer. This is 
particularly true in an interfacial region or near a free sur- 



face, where material properties can change significantly 
over short distances. 

In this work, we are interested in the local elastic 
constants and surface melting of thin crystalline films. 
Specifically, we present a formalism in which the local 
elastic constants are calculated at precise planes in the 
system, as opposed to small volumes or slabs. In the 
bulk, the calculated local elastic constants are verified 
by averaging over the entire system and comparing the 
results to the bulk value. The local elastic constants in 
the film are verified by comparing the local stress pro- 
files that arise from uniaxial strain and those calculated 
directly from the elastic constants. 



II. THEORY 

In a homogeneous material, applying a homogeneous 
strain necessarily results in a homogeneous stress. The 
stress is given by 



(1) 



where Cijki is the bulk elasticity tensor, and where the 
indices represent the cartesian coordinates in three di- 
mensions. 

When a homogeneous strain is applied to an inhomoge- 
neous system, the resulting stress is also inhomogeneous. 
The local stress is then given by 



o-y(r) = Cy7ci(r)e/fc, 



(2) 



where C'ijki (r) is the local elasticity tensor. The relation- 
ship between the local and bulk elasticity tensors can be 
written as 



(3) 



where V is the volume of the system. 
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The bulk elasticity tensor can be expressed in terms of 
the fluctuations of stress according to M 



V 
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where p is the density, 6ij is the Kronecker delta, is 
the pressure tensor and Bijki is the so-called Born term. 
The pressure tensor is given by 
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The potential energy between interaction sites a and b 
is denoted by Uab, Tab is the distance between them, pa 
and nia are the momentum and mass of site a respec- 
tively, and the prime indicates a derivative with respect 
to Tab- The Born term is related to the first and second 
derivatives of the potential energy of interaction by 
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In this work we focus on the elastic properties of thin 
films having a planar symmetry; the films are assumed 
inhomogeneous only in the direction perpendicular to the 
film, i.e. z. Equation (Q) must therefore be modified to 
calculate the elasticity tensor at precise planes within the 
system, without need for bins or small volumes. To this 
end, we use the method of planes (MOP) Q and obtain 
an expression for the local elasticity tensor. 

The first term in Eqn. (Q) is the ideal gas contribution 
to the elasticity tensor. The kinetic energy is homoge- 
neously distributed, even in inhomogeneous systems, and 
the temperature is independent of z. However, the den- 
sity can vary in the z direction. The density profile, p{z), 
could be calculated by dividing the system into many 
small bins and counting the average number of particles 
per unit volume in the bins. The density would then ex- 
plicitly depend on the size of the bins used. Alternatively, 
one can use the fact that for a free standing film the to- 
tal normal pressure, Pzz = p{z)kBT + P^^{z), is constant 
throughout the system pO| . In vacuum, we then have for 
the density profile 



(7) 



where P^zi^) is the configurational contribution to the 
local pressure tensor. The local pressure tensor is the sum 
of ideal and configurational terms, and can be calculated 
according to |lOl 0| 
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where A is the cross-sectional area of the film and Q is 
the Heaviside step function. The first term in Eqn. (^) 
can then be written for inhomogeneous systems as 



2p{z)kBT [SiiSjk + SikSji] 



(9) 



The second term in Eqn. arises from bulk stress 
fluctuations; it accounts for the non-zero temperature 
contribution to the elastic constants. We are interested 
in relating the local stress, (t(z), to a bulk homogeneous 
strain. Therefore, instead of including the bulk-stress 
bulk-stress correlation, we use the correlation between 
the local-stress and the average bulk stress. The second 
term can then be written as 
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Note that the volume V in Eqn. (|1C|) cancels that in 
Eqn. ^ and there also is no explicit volume dependence 
in the MOP expression for Pij{z). 

The Born term, Eqn. (|^), can be calculated at planes 
using the MOP in the same way the local stress is de- 
termined, i.e. Eqn. (||). We have for the Born term in 
inhomogeneous systems 
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As before, this expression does not depend on the volume 
of the system or the (arbitrary) size of a bin. 

The flnal expression for the local elasticity tensor in 
inhomogeneous systems with planar symmetry is given 
by 



(12) 



We emphasize again that this expression for Cijki (z) is 
valid for inhomogeneous systems and is an average only 
over a cross section (a plane) of the system, not a discrete 
volume. It therefore relates the local stress (at z) to a 
homogeneous strain. Also note that Lutsko et al. p| have 
presented a derivation for the local elasticity tensor, but 
they averaged over a sub- volume in order to facilitate the 
computations. It can be seen that by integrating over 
the entire system, one recovers the bulk elasticity tensor, 
Eqn. (^). We also note that this expression does not 
require the use of any dynamic variables but only requires 
ensemble averages taken from system configurations. It 
therefore is useful in either molecular dynamics or Monte 
Carlo (MC) simulations. 

We note that other valid definitions of the local stress 
tensor have been presented 
extensively in the literature 



12, 14 1 and discussed 
L5|. These definitions 
would in principle lead to different expressions for the 
local elasticity tensor. Regardless of the definition, one 
should expect to recover the bulk elasticity tensor after 
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averaging over the entire system. The definition used 
in this work is that of Irving and Kirkwood . This 
definition was chosen here because it has been shown to 
be a physically valid stress tensor and it can be used 
in MC simulations. 



III. SIMULATIONS 

To demonstrate the calculation of the local elastic- 
ity tensor, we employ the widely used nearest-neighbor 
Lcnnard-Jones (NNLJ) fee crystal model |l^, 

pl| . In what follows, all results will be reported in 
dimensionless Lennard-Jones units. 

A bulk system consisting of 1000 particles with peri- 
odic boundary conditions in all three dimensions was in- 
vestigated first. This system was simulated in the NVT 
ensemble at a temperature of T = 0.3 using a simple MC 
method. The density was chosen such that the average 
bulk pressure is zero. The center of mass of one atomic 
layer was held fixed at 2 = 0. The average bulk elastic 
constants for this system have been calculated previously 
and are listed in Table |[ Reported elastic constants, 
stresses and strains are represented using Voigt nota- 
tion For bulk fee systems, there are three groups 
of non-zero, independent elements of the elastic constant 
matrix 
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Second, we also consider a free standing film of 450 
particles. The free surfaces correspond to the (001) face 
of the fee crystal. This system was also simulated in 
the canonical ensemble using a conventional MC method. 
In this case the cross sectional area was held constant 
with the same dimensions as the bulk system. The film 
had nine atomic layers parallel to the free surfaces. The 
temperature was the same as in the bulk, i.e. T — 0.3. 
The center of mass of the film was held fixed at z = 0. 
For an fee film with free surfaces normal to the z-axis, 
there are six groups of non-zero, independent elements of 
the elastic constant matrix 



(14) 



Local properties of the film and the bulk system were 
calculated from Eqn. ( p^ at planes of constant z, with 
each plane being separated by a distance of 0.02 in both 
the thin film and bulk systems. The average elastic con- 
stants were also calculated in the bulk system from Eqn. 
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FIG. 1: Cii(z) (solid line) as a function of z for the bulk 
system from Eqn. (p^. The density profile, pb{z), is shown 
as the dotted line. 
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FIG. 2: 6*21(2) (solid line) as a function of z for the bulk 
system from Eqn. (|l2|). The density profile, pb{z), is shown 
as the dotted line. 



An additional simulation of the thin film was per- 
formed in which a homogeneous, tensile, uniaxial strain 
was applied in the x-direction, ei = 0.01005. The strain 
is defined as 1231 



(15) 



where is the length of the simulation cell in the x- 
dircction, and L*J is its original length. Since the strain 
is homogeneous, it is known that the average strain in 
a plane of atoms parallel to the free surface is equal to 
the applied strain [Q. The resulting stress profiles were 
then calculated using Eqn. (||). The stress profiles were 
also calculated directly from the elastic constants using 
Eqn. (I). 



IV. RESULTS 

The local elastic constant profiles for Cii{z) and C2i{z) 
in the bulk system are shown in Fig. |l] and Fig. ^ re- 
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Ref [25] 


Eqn. (4 


Cii 


43.35 


43.22 


C21 


19.01 


19.45 


C44 


22.50 


22.60 



JC-j{z)dz 



43.37 
19.22 
22.35 



TABLE L Values of the three independent elastic constants 
of the bulk fee crystal in dimensionless Lennard- Jones units. 
The last column is the average value of eight atomic layers of 
the (001) face in the bulk from Eqn. (|l|). 
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FIG. 3: C(-i{z) (solid line) as a function of z for the film 
system from Eqn. (p^. The density profile, Pf{z), is shown 
as the dotted line. 



spectively. The density profile, pb{z), is also shown in 
these figures. Each peak in the profile at C^i « 175 
and C21 « 150 corresponds to the center of mass of each 
atomic layer. Each minimum at « 22.5 and C21 ~ 
corresponds to the midpoint between each atomic layer. 
The local elastic constant profile for C44 is similar to C21 
and is not shown. 

The average bulk elastic constants can be calculated 
from the local constants using 



(16) 



where Z is the width of the system. In this example, 
we set Z to the width of eight atomic layers in the bulk. 
The average bulk elastic constants calculated from Eqn. 
( p^ ) are given in Table ^. The bulk elastic constants from 
Eqn. (0) and the literature values are also given in 
Table |7 All three values for each elastic constant agree 
well with one another. 

The local elastic constant profiles for C{-^ (z) and C/i {z) 
in the thin film are shown in Fig. ^ and Fig. ^, respec- 
tively. The density profile, Pf{z), is also shown in these 
figures. The peaks corresponding to the atomic layers 
in the center of the film (z = 0) have approximately 
the same maximum values as in the bulk system, i.e. 
C/^ (0) « C^(0). However, the minimum values between 
each layer near the center of the film are less than in the 
bulk. Interestingly, Cli{z) exhibits negative values be- 
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FIG. 4: Ci^iz) {solid line) as a function of z for the thin film 
from Eqn. (p^. The density profile, pf(z), is shown as the 
dotted line. 
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FIG. 5: C(i(z) as a function of 2 for a film with 17 atomic 
layers. 



tween each layer. The meaning of these negative values 
is discussed below. 

The peak values of the elastic constants decrease from 
the center of the film as the free surfaces are approached. 
The profiles also become broader near the free surfaces. 
The decrease of the local elastic constants is an indication 
of the enhanced atomic mobility at the surfaces. 

In order to investigate the effect of film thickness, a 
film consisting of 17 atomic layers was also simulated. 
Figure ^ shows the local elastic constant profile for C(i (z) 
in the film with 17 layers. The effect of the free surface is 
limited to the first two atomic layers for both this system 
and that shown in Fig. ^ The elastic constant profiles 
for both film thicknesses are consistent with one another. 

The local stress profiles are shown in Fig. || for cri(z) 
and in Fig. ^ for cr2(z) after a homogeneous uniaxial 
strain was applied. The local stress profiles in the film 
were calculated from Eqn. (|^) using the local elasticity 
tensor measured at zero strain. The local stress profiles 
were also calculated in the strained film as 



(17) 
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crystal. The interaction is described by the truncated LJ 
potential given by 



FIG. 6: Profiles for cti in the thin film after a small homoge- 
nous uniaxial strain, ei , is applied. The solid line is calculated 
from the elastic constants and Eqn. and the dotted line 
is calculated directly from the simulation using Eqn. (|q). 




FIG. 7: Profiles for (T2 in the thin film after a small homoge- 
nous uniaxial strain, ei , is applied. The solid line is calculated 
from the elastic constants and Eqn. and the dotted line 
is calculated directly from the simulation using Eqn. m). 



where we used Eqn. (g) for Pi{z). The results are shown 
in Fig. H and Fig. 0. The fact that the two methods for 
calculating the local stress profiles give the same result 
is reassuring and demonstrates that the response to the 
applied strain is linear. For a2{z), we find that the tensile 
(positive) uniaxial strain in a;-direction causes a negative 
stress in the region between the atomic layers (Fig. |^). 
This is directly related to the negative elastic constants 
seen in C|j(z). 



V. SURFACE MELTING 

In order to study the melting behavior of a thin crys- 
talline film, we adopt the model for argon used by Eerden 
et al. p^. As before, we study the (001) surface of the 



Uab =4.569e 



X exp 



0.25g 
rab - 2.5cr 



(18) 



Eerden et al. report the bulk elastic properties for this 
system. 

Consistent with Eerden et al., we have 32 atoms in 
each layer of the crystal and use films consisting of 16 
layers. The elastic constants are calculated at planes in 
the top half of the film separated by a distance of dz — 
0.02. At each temperature, the lateral dimensions of the 
simulation cell were taken from the average size of a bulk 
simulation cell at zero pressure. The surface of the film 
was aligned perpendicular to the z-axis and the center of 
mass was fixed at 2; = 0. 

Analogous to the definition of the average lateral shear 
modulus for a slab between z and z' {^"^[z, z'\) by Eerden 
et al., we define the local lateral shear modulus at a plane 
z as 



/^r(2) =^ XI 51 ['^"/3/5"(^) 

a.—x,y (3—x,y 



(19) 



Note that this definition is a projected (onto the xy- 
plane) version of the usual shear modulus for an isotropic 
solid . Since we are interested in the melting behavior 
of an anisotropic solid, another useful definition of the 
local lateral shear modulus is 



<(z)=C66(z). 



(20) 



As the crystal nears its melting point it becomes less 
anisotropic and we expect /x^ to approach /if. The melt- 
ing point is defined here as the temperature at which /i^ 
and /zf vanish. 

The shear moduli as a function of position in the film 
are shown in Fig. |^ at four different temperatures. The 
density profiles at these temperatures are shown in Fig. ||. 
The density profile has units of and its integral over 
the entire system, A p{z)dz, gives the total number 
of particles in the film. In the following discussion, we 
will refer to the layers starting with the surface layer as 
layer- 1, layer-2, etc. 

The behavior at the lowest temperature (Fig. ||a.), 
T = 0.4, is similar to that of the NNLJ film, having bulk 
behavior in the center of the film and decreasing moduli 
in the layers near the surface. The difference between /if 
and /if reflects the fact that the crystal is anisotropic, 
even in the layer nearest to the surface, layer- 1. This is 
also evident in the density profile (Fig. ga.) where all the 
atomic layers of the crystal are separated by regions of 
empty space. At temperatures above T = 0.4, isolated 
atoms have sufficient energy to escape layer- 1 and occupy 
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FIG. 8: The shear moduli profiles in a thin film of argon at 
four temperatures. The solid line is the fi^ and the dotted 
line is fif . 




FIG. 10: The average lateral shear moduli in a thin argon 
film as a function of temperature. 
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FIG. 9: The density profiles in a thin film of argon at four 
temperatures. 

positions outside the film (layer-0) . This additional layer 
can be seen as the very small peak centered at z = 6.82. 
The additional layer, however, has a zero shear modulus. 

At r = 0.5, the shear moduli (Fig. of each layer 
have decreased from those at T — 0.4, indicating a soft- 
ening of the crystal. The difference between fj,f and /i^ 
has decreased considerably in layer-1, indicating nearly 
isotropic behavior near the melting temperature. The 
number of atoms which escape from layer-1 has increased, 
indicated by the larger peak or shoulder in the density 
profile (Fig. |b.) at z = 6.82. 

In Fig. |c. and Fig. |c. (T = 0.575), the behavior 
of the layers near the surface has changed significantly. 
Both /if and /i^ are essentially zero (indicating isotropy 
and melting) at layer-1 even though the density profile 
exhibits some remaining structure in that region. Be- 
tween layer-1 and layer-2 and between layer-2 and layer- 
3, the density is non-zero yet the shear modulus is zero. 
A small amount of argon exists as a fluid between these 
layers. Layer-0 contains even more atoms at this tem- 



perature and shows a flat density shoulder which decays 
to zero, indicating a loss of structure at the fllm- vacuum 
interface. At temperatures just below T = 0.575 and 
above, the delineation of layer thickness becomes am- 
biguous and the use of a layer-averaged shear modulus 
becomes questionable. The method of planes proposed 
here eliminates that ambiguity. 

At T = 0.6, the shear modulus of the entire film is 
zero and the density profile is flat. The film is a liq- 
uid throughout and has none of the structure originally 
present in the crystalline film at lower temperatures. 

An average shear modulus for the surface layer can be 
defined by integrating the profiles in Fig. |^. The average 
shear modulus is given by 

r = ""°%"(^)d^. (21) 

A layer thickness Az must first be defined in order to 
perform the integration. We arbitrarily choose Az for the 
surface layer to be the distance between the peaks in the 
density profile of layer-1 and layer-2 at each temperature. 
In Eqn. (pl|), 2;,„i„ is the location of the minimum density 
between layer-1 and layer-2 and Zmax = -Zmin + Az. 

The results for Ji? and Ji^ of the surface layer are shown 
in Fig. |l^ for temperatures up to T = 0.575. Both 
shear moduli decrease sharply with increasing temper- 
ature and vanish at T = 0.575, indicating melting of the 
surface layer. This is in agreement with the literature 
value of r = 0.576. The bulk melting temperature is 
Tb = 0.601 |§. 

VI. CONCLUSION 

A formalism for calculation of the local elastic con- 
stants in inhomogeneous systems based on the method 
of planes has been presented. Unlike previous work, this 
formalism does not require the partitioning of the system 
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into a set of finite volumes or "slabs" over which average 
elastic constants are calculated. As an example appli- 
cation of the technique, Monte Carlo simulations of the 
nearest-neighbor Lennard-Jones fee crystal in the bulk 
and in thin film geometries have been presented. 

The local atomic structure of the crystals was evident 
in the local elastic constants calculated at precise planes. 
In the thin film, the elastic constants are decreased from 
the corresponding bulk values, especially near the free 
surfaces. This decrease near a free surface is expected 



to give rise to apparent deviations from bulk continuum 
behavior in thin films and nanoscopic structures. 

The melting behavior of argon in a thin film was also 
investigated within the context of this formalism. Results 
show how the shear modulus profile of the surface layer 
of atoms vanishes below the melting temperature of the 
core of the film. Below the melting temperature of the 
film, the free surface allows sufficient thermal motion for 
the surface atoms to reach an isotropic liquid state prior 
to the bulk of the film. 
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